Non-classical velocity statistics in a turbulent atomic 
Bose— Einstein condensate 
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^ Abstract 
< 

In a recent experiment Paoletti et al. (Phys. Rev. Lett. 101, 154501 (2008)) monitored 
the motion of tracer particles in turbulent superfluid helium and inferred that the velocity 
components do not obey the Gaussian statistics observed in ordinary turbulence. Motivated by 
their experiment, we create a small turbulent state in an atomic Bose-Einstein condensate, which 
^ enables us to compute directly the velocity field, and we find similar non-classical power-law tails. 

Our result thus suggests that non-Gaussian turbulent velocity statistics describe a fundamental 
^ property of quantum fluids. We also track the decay of the vortex tangle in the presence of the 

1 thermal cloud. 
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Quantum turbulence (a dynamic tangle of discrete, reconnecting vortices) has been 
studied in superfluid "^He P], superfluid ^He-B [2], and, more recently, in atomic Bose- 
Einstein condensates (BEC's) |3l IH E]- The defining property of these quantum fluids is 
that the superfluid velocity fleld is proportional to the gradient of the phase of a complex 
order parameter, so the circulation is quantized in unit of the quantum k (the ratio of Planck 
constant to the mass of the relevant boson). Therefore, whereas in classical ordinary fluids 
(such as air or water) the rotational motion is unconstrained, the velocity fleld around a 
quantum vortex decreases strictly as /t/(27rr) where r is the distance away from the vortex 
axis. 

Recent work has revealed that there are remarkable similarities between classical 
turbulence and quantum turbulence: the same Kolmogorov energy spectrum in continuously 
excited turbulence [6], the same temporal decay of the vorticity in decaying turbulence [7], 
the same pressure drops in pipe and channel flows [8], and the same drag crisis behind a 
sphere moving at high velocity [9J. These results have attracted attention to the problem 
of the relation between a classical fluid (which obeys the Euler or Navier-Stokes equations) 
and a quantum fluid [lO] , stimulating apparently simple questions such as whether classical 
behaviour is merely the consequence of many quanta. A recent experiment by Paoletti 
et al. [11] in superfluid ^He has shed light onto this problem. Using a new visualization 
technique based on solid hydrogen tracers, Paoletti et al. found that the components of the 
turbulent velocity field do not obey the usual Gaussian distribution which is observed in 
ordinary turbulence [TB], but rather follow power-law like behaviour. 

The main aim of this Letter is to answer the question of whether non-Gaussian velocity 
statistics is a fundamental property of turbulence in a quantum fluid. To achieve this aim, 
we move from ^He to an atomic BEG. The agreement that we flnd between these two 
distinct systems means that power-law behaviour is indeed typical of quantum fluids, in 
stark contrast with classical turbulence. 

To create a vortex tangle in a harmonically conflned atomic BEG we choose the technique 
of phase imprinting [13]. Although phase imprinting is not the only method which can be 
used to generate turbulence [3 [12], it is speciflc to BEG's, and can in principle be exploited 
to non-destructively generate tangles engineered to be isotropic or polarized; this gives the 
technique its own numerical advantages over others. Having generated a vortex tangle, we 
also study its decay at zero temperature by evolving the three-dimensional (3D) Gross- 
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Pitaevskii equation for realistic experimental parameters, and also discuss the effect of the 
thermal cloud on this decay. 

An ultra-cold dilute atomic gas confined in a spherical trap is described by the Gross- 
Pitaevskii equation (GPE) for the complex wavefunction \1'. It is convenient to consider the 
GPE in dimensionless form as 

9 ^ 1 2 , 1 2 , ^, ,.,2 



where r is the position vector; tjj = a^^'^N~^^'^'^ is dimensionless and satisfies the constraint 



/ dr lipl = 1. In writing Eq. (1 ) we use the harmonic oscillator length Oo = yh/ (muj) as the 
unit of distance and the inverse trapping frequency l/tu as the unit of time, where m is the 
mass of one atom, C = AnNa/ao is the dimensionless measure of the interaction between 
bosons, is the number of atoms and a is the scattering length. 

Finite temperature effects can be simulated by replacing i at the L.H.S. of Eq. ([!]) with 
(i— 7), where 7 models the dissipation arising from the thermal cloud. Although such a model 
can be justified from first principles [H], one often resorts to the use of a (phenomenological) 
constant dissipation term [T5] (as done here), as it is much less computationally intensive 
while still giving a qualitatively accurate prediction of the relevant physics. 

We employ realistic experimental parameters for a ^^Na condensate (a = 2.75 nm) of A^ = 
10^ atoms and trapping frequency u = 271 x 150 Hz, giving C = 2.019 x 10'^, = 1.71 /im 
and 1.06 ms as time unit. In a homogeneous condensate the healing length ^ is estimated 
by balancing the kinetic energy per particle and the interaction strength, which implies the 
(dimensionless) ^ = (2Cn)~^/^ where n = l^/'p is the condensate density. In a harmonic trap, 
n is position dependent, so we estimate the average ^ using the mean density; for example 
at t = to (Fig. 0, (^) = 1-35 x 10"=^ and = 0.43. 

Eq. ([l| is evolved pseudo-spectrally via XMDS [I7j in 3D using a 4*^^ order Runge-Kutta 
method and periodic boundary conditions. The spatial domain |x|,|?/|,|2;|<8is discretised 
on a 128 x 128 x 128 grid with timestep At = 10~^. We have verified that the results are 
independent of spatial and temporal stepsize. 

The initial condition is created by phase-imprinting a grid of 16 vortices oriented parallel 
to the Cartesian axes in a staggered way (no vortices intersect), and propagating this 
configuration for a short period in imaginary time (which is equivalent to substituting t 
with —it in Eq. ([T|) while continuously re- normalizing ip, until the density adjusts to reveal 



-5 -5 



(a) 




°0 2 4 6 8 10 

(b) 



FIG. 1: (Color online) (a): Generated turbulent state within the condensate edge (blue shading), 
at to = 1-9 and 7 = 0; the maximum density (used to determine the condensate edge) is n = 
4.388 X 10^'^. (b) Evolution of normalized total length L{t — t())/LQ for 7 = (black asterixes, *) 
0.015 (purple pluses, +), 0.03 (blue circles, o) and 0.06 (red diamonds) corresponding to the same 
initial phase imprinting on a 128'^ grid; the inset also shows the initial increase of L for t < Iq and 
the exponential fit for 7 = (solid green line). 

the desired vortex structure. This is then propagated in real time. In the initial evolution 
period, the vortices reconnect and become excited, resulting in an increase in the total vortex 
line length L. This increase continues up to some time to, when the L{t) achieves its peak 
value Lq, and a maximally tangled turbulent state has been generated (Fig. [l|^a)). At this 
point the dimensionless vortex line density (vortex length divided by volume) is 0.79. The 
vortex line length evolution for t > to is shown in Fig. [l](b). Since no energy is injected 
into the system the turbulence is not sustained but decays over a time scale of the order 
t ~ 10, as vortices further reconnect, break up and decay into sound waves [18], or leave 
the condensate. Our calculation lasts longer (up to t = 27) at which time one long vortex 
is left at the centre of the condensate and some shorter vortices, roughly aligned in the 
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same direction, surround it closer to the condensate's edge. Large volume oscillations of 
the condensate at 1/3 of the trap frequency are observed. Although the amplitude of these 
oscillations decreases as the turbulence decays, a method of exciting turbulence without 
inducing volume oscillations would be preferable for experiments. 

At each time t we attribute a vortex core length to every point in the condensate at which 
the real and imaginary parts of ip crosses zero (therefore defining the vortex axis) [19]. The 
total vortex length L{t) is the sum of all identified vortex points within the condensate 
edge, which, throughout this work, is defined as the outermost points at which n drops 
below 25% of the maximum density. We find the decay of the vortex length is better fit by 
an exponential of the form L{t — to)/Lo ~ exp(— c(t — to)) (solid line in Fig. 1(b)), rather 
than 1/L{t — to) (x t — to found in p[]. Fitting over 10 time units, yields consistent values 
for different grid sizes (c = 0.151 ± 0.005 and 0.144 ± 0.007 for gridsizes of 128^ and 256^ 
respectively, with corresponding line lengths Lq = 406 and 378). We have also checked that 
the decay of L{t) is largely insensitive to the initial imprinted vortex configuration obtained 
by imprinting extra vortices. 

To obtain a qualitative understanding of the effect of the thermal cloud, we have repeated 
the above calculation with dissipation parameters 7 = 0.015, 0.03 and 0.06. As anticipated, 
the role of temperature is to induce faster decay of the turbulent state, leading respectively 
to c = 0.252 ± 0.007, 0.274 ± 0.006 and 0.340 ± 0.018 (with corresponding Lq = 369, 339, 
and 300, due to the damping in the initial period t < tg)- 

Our next step is to analyse the turbulent (dimensionless) superfiuid velocity field, which 
we compute directly from the definition v(r) = [tp* \/ ip — ip \/ / {2i\ilj\'^). (the derivatives 
of ip being obtained spectrally). The calculation is tested against the expected azimuthal 
velocity profile of a single vortex set along the z axis, which is Vg = K/(27rr) where k = 27r 



is the quantum of circulation. Using the Madelung representation ip = ^nir) exp[iip{r)], 
yields v(r) = Vip{r), showing that the velocity depends only on the condensate phase (p. 

We calculate the probability density function (PDF) of each Cartesian velocity component 
Vi {i = X, y, z) for 7 = and compare it to the form that a Gaussian PDF (gPDF) of the 
velocity would take: 



where cr, cr^ and /i are the standard deviation, variance and mean of the velocity distribution 
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FIG. 2: (Color online) Left: 3D log-linear velocity PDFs corresponding to Fig. 1 (7 = 0): 
logioPDF(fj)] vs. velocity component Vi, for Vx (blue circles, o), Vy (red triangles, v) and Vz 
(green asterixes, ) , yielding following power law coefficients respectively b = — 3.27ib0.04, — 3.54ib 
0.06, —3.57 ± 0.06. Corresponding logxo[gPDF(w.j)] plots shown by black dotted (• • •, Vx), dash- 
dotted (-., Vy) and uniform v^) lines, which almost overlap. The velocity components are only 
sampled within the condensate edge (defined at 25% of the maximum density) thus excluding the 
outer region where the density drops to zero. Corresponding statistical values: af = 4.1, 3.6, 4.2 
and fli = (4.6, 4.5, 4.6) x 10~^. Right: Log-log plot of the same PDFs; the solid line is the power 
law fit to Vx- 

(see Fig. |2|. In evaluating the gPDF, we use the mean, standard deviation and variance 
of the actual velocity PDF. Power-law dependence of PDF(fj)oc is found from log(PDF) 
vs log(i;j) plots of the positive velocity components, with —3.57 < b < —3.27 in all three 
directions. 

It is apparent that our velocity statistics are non-Gaussian, consistently with the high 
velocity tails found experimentally by Paoletti et al. [Hj in turbulent superfiuid ^He (their 
h is —3, slightly less than ours). It is well-known that the PDF's of velocity components 
in ordinary turbulence are Gaussian [16]; our result thus supports the idea (put forward in 
Ref. [TT]) that non-Gaussian velocity statistics are the unique signature of the quantized 
nature of the vortices. 

We find that deviations from Gaussian behaviour are less pronounced (but still noticeable) 
if we omit sampling the velocity very near the axis when the density is less than a prescribed 
cutoff. Unlike a classical vortex, there are not individually distinguishable atoms which spin 
about the vortex axis, and the velocity field is entirely defined by the macroscopic phase 
(y9(r) irrespective of the density Ti(r), so our procedure is justified. The small wiggles in the 
PDF for values close to zero may be a measure of the anisotropy of the vortex tangle in our 
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condensate and would not feature in velocity PDFs of a truly isotropic state of quantum 
turbulence in a larger condensate, or they could be due to the large volume oscillations of 
the 3D condensate. 

The observed non-Gaussianity of the velocity statistics holds during the decay of the 
vortex tangle, suggesting that it is not necessarily caused by vortex reconnections, whose 
frequency depends on the vortex line density |20]. To verify that our result is general and 
does not depend on vortex line density or reconnections, we repeat the calculation in a two- 
dimensional (2D) BEC; such a condensate is created when the axial trapping frequency, uOz, 
is much greater than the radial trapping frequency, cj,., freezing out motion along z. The 
condensate in the radial plane is also described by Eq. with C — > C2d = 2^/271 aN / az-, and 
4'2d = arN^^^'^"^ which relates the dimensionless wavefunction to the dimensional 2D 
condensate wavefunction, \l/2d, and = \jTi/{muJr) (here = j [raoj are the harmonic 
oscillator lengths in the radial and z directions). To induce turbulence in a 2D BEC we 
imprint equal numbers of oppositely charged vortices randomly located and aligned along 
the axial direction [3]. Here we simulate a 2D ^^Na Bose gas with = 10'' atoms and 
condensate aspect ratio uJz/oJr = 20 (with ujr = 7.5 x 2tt Hz), implying C2d = 8.06 x 10^. 
The spatial computational domain |?/| < 25 is discretized on a 1024 x 1024 grid with 
timestep At = 10~^. Choosing our initial monitoring time arbitrarily at to = 4.4, we 
have 86 vortices in the dimensionless condensate area 752, corresponding to a dimensionless 
vortex line density (number of vortices per unit area) of 0.11. We find a mean radial 
density {n2d) = 1-23 x 10^"^, and an estimated healing length ^2d = {2C2dn2d)~^^'^ = 0.225. 
Proceeding as in 3D, the velocity PDF's and gPDF's of the 2D condensate are shown in 
Fig. 3 for 7 = 0; we find again non-classical, non-Gaussian velocity statistics, with high 
velocity tails. This PDF lacks the small wiggles close to zero velocity observed in the 3D 
case, as the method of phase imprinting 42 positive and 42 negative vortices along the z axis 
in randomly selected locations across a much larger condensate generates a more isotropic 
turbulent state with negligible oscillations of the total radial area. As in 3D, we find a 
power-law correlation with velocity, PDF(fj) oc with —3.18 ± 0.08 for positive Vx and Vy 
velocity components respectively. Unlike the 3D case, vortex reconnections have played a 
negligible role in the state which is measured because during the measured timescale the 
number of vortices remains approximately the same as initially imprinted. 

In conclusion, we have shown (under realistic experimental conditions) how phase 
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FIG. 3: Density (left) and velocity PDFs (right) for a 2D BEG at to = 4.4 and 7 = 0. The 
maximum density is n = 2.057 x 10~^. Left: 86 vortices (core radius ~ 2.66) can be identified. 
Right: Plots of logioPl-'F(fj)] and the corresponding log;^o[gPI^F(?;i)] vs where the Vx^ Vy 
velocity components (sampled within the condensate edge) are given by blue circles (o) and red 
triangles (y); corresponding gPDF's displayed by black dotted (• • •) and dash-dotted (-.) lines. 
Statistical values: = ay = 2.0, jlx = 5.8 x 10^^, fiy = 1.3 x 10^^. 

imprinting of a staggered array of straight non-intersecting vortices in a harmonically trapped 
ultra-cold atomic gas can be used to generate and study quantum turbulence in atomic 
EEC's. Finding velocity statistics similar to ^He in these relatively small systems means 
that atomic EEC's can be used in the study of turbulence. We have also found that the 
decay of turbulence is faster in the presence of a thermal cloud, which we have modelled in 
a simple way. 

Our main result is that the statistics of the turbulent superfluid velocity components are 
non-Gaussian, and exhibit power law tails similar to what has been recently found in an 
experiment with turbulent superfluid ^He [TT]. We confirmed our 3D result by repeating the 
calculation in a 2D condensate with a greater number of vortices, but lower dimensionless 
vortex line density, in a system in which vortex reconnections played virtually no role. 
Despite the huge difference in the ratio of intervortex spacing to vortex core radius in atomic 
EEC's (~ 3 here) and in ^He (10^ to 10^), we observed the same non-Gaussianity of the 
velocity, in contrasts with the known Gaussianity of the velocity in classical turbulence. Our 
result thus suggest that power-law tails are a general feature of the constrained 1/r velocity 
fields which characterizes a quantum fluid. 

The statistics of velocity components are thus macroscopic observables which distinguish 
between classical and quantum turbulence. It is worth remarking [21] that another such 
observable is the pressure spectrum, which, in 3D quantum turbulence, due to the 1/r 
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velocity field and the singular nature of the vorticity, should obey a law (where k is 
the magnitude of the wavevector) rather than the k~'^^^ scaling which corresponds to the 
classical Kolmogorov k~^^^ spectrum of the energy. 

This research was supported by EPSRC grant EP/D040892/1 and by the Merit Allocation 
Scheme of the National Facility of the Australian Partnership for Advanced Computing. 
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